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Abstract 

For sampling multiple pathways in a rugged energy landscape, we propose a novel action-based 
(T) ■ path sampling method using the Onsager-Machlup action functional. Inspired by the Fourier-path 

o. 

integral simulation of a quantum mechanical system, a path in Cartesian space is transformed into 
that in Fourier space, and an overdamped Langevin equation is derived for the Fourier components 



to achieve a canonical ensemble of the path at a finite temperature. To avoid "path trapping" 
around an initially guessed path, the path sampling method is further combined with a powerful 
sampling technique, the replica exchange method. The principle and algorithm of our method is 
numerically demonstrated for a model two-dimensional system with a bifurcated potential land- 
scape. The results are compared with those of conventional transition path sampling and the 
equilibrium theory, and the error due to path discretization is also discussed. 
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I. INTRODUCTION 



To understand the functions of biomolecules (chemical reactions, conformational change, 
Ugand-binding, protein-protein aviation etc.), detennhnng reaction paths is an essentia, 
step However, conventional molecular dynamics (MD) simulations often fail to find 
appropriate reaction paths even with huge computational resources j^j. This is because 
such reactions cannot be simulated within computationally feasible times, and the sampling 
efficiency of the reaction can be low especially for large systems. Handling such a rare event 
is a common issue not only for biomolecular systems such as proteins but also for many 
physical or chemical systems. 

So far, many methodologies have been devised for path search or path sampling in 
biomolecular systems The most primitive but powerful method may be umbrella sam- 
pling after defining certain reaction coordinates such as the radius of gyration or the number 
of native contacts 

a a. 

Alternatively, metadynamics has been used to compute free en- 
ergy surfaces with several reaction coordinates [5] . However, all these methods only work for 
molecular processes characterized by limited numbers of reaction coordinates that are known 
a priori, and one often fails to describe complex and highly multi-dimensional processes in 
biomolecular systems. As the computational effort grows exponentially with the number 
of reaction coordinates (dimensions) increases, it is hopeless to work on the full detail of a 
high- dimensional free energy landscape. Steered or targeted molecular dynamics is another 



powerful and well-established path search method [6|, |7[ • In combination with the Jarzynski 
equality UQ, a free energy difference can be calculated by steered MD. However, since 
the reaction path obtained only follows a predefined direction, it may sometimes result in 
producing paths with unphysically large deformations of molecules. 

More direct approaches to deal with complex multi-dimensional reactions are based on 
minimization of a physically defined action functional associated with a reaction path. El- 
ber and Karplus JjJ proposed an innovative approach to use a line integral along a path 
with some constraints, in search of a minimum energy path such as the intrinsic reaction 
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path 
Other re 
method 



14{ . This method was successfully applied to peptides [15| and water systems 



16] . 



ated methods searching for a minimum energy path include the nudged elastic band 
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string method 



the conjugate peak refinement method [19|, [20(, and the zero-temperature 
211 ] . These methods can be used to locate transition state geometries (saddle 
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points of the potential energy surface, at which barrier crossing can occur with a minimum 
amount of activation energy), but they do not contain any dynamical information by them- 
selves, since the transition state geometry is merely a static feature of the potential energy 
landscape. Some kinetic information is retrieved with the use of the discrete path sampling 
method [22| under the assumptions of transition state theory 23f| , but the effect of thermal 
fluctuations, which should be particularly important for biomolecular systems, is not com- 
pletely captured by these approaches. To include dynamical and thermal effects, transition 
path sampling (TPS) has been proposed and developed by Chandler and coworkers 24. 125]. 



The basic idea of TPS is that a system is described bvan ensemble of paths instead of a 
single path. (For a similar and heuristic approach, see [261] .) When an overdamped Langevin 
equation is assumed to describe the dynamics, the weight for a path can be written as 



P[x(t)] oc e 



-S[x(t)]/2D 



(1) 



where 



S[x(t)} 



AU 1 
- + 2 



ds 



.o i ldliy 2Dd 2 U 



C dx J 



is called the Onsager-Machlup (OM) action "functional" 
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(2) 



pedagogic discussions, see, for example, 



C dx 2 

28j (for its derivation and 



29] ) because it contains all the information of 



the whole history of a path x(t). Here, AU = U(xsn) — U(xi n i) is the potential energy 
difference between an initial state with x ini = x(0) and a final state with x^ n — x(t), D is the 
diffusion constant, which is related to a friction constant ( and the absolute temperature T by 
Einstein's relation D = ksT '/( (ks is the Boltzmann constant). TPS intends to sample many 
paths according to the statistical weight above. Combining conventional MD simulations 
and Monte Carlo (MC) moves, shooting and shifting [24] . TPS has been successfully used 
for rare and fast chemical reactions. 

However, in biomolecular reactions, such as folding or conformational changes, the pro- 
cess is often not only rare but also slow and diffusive. Protein folding can occur at least 
more than microsecond timescales, which is still beyond the ability of current "brute-force" 
MD simulations using "standard" computers and algorithms. This is also the case for con- 
formational change of large protein systems with several domains, which is relevant to the 
understanding of the functions of proteins [1]. Therefore, alternative approaches to sample 
path ensembles should be established, and novel hardware strategy {2]] and novel theoretical 
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algorithms including transition interface sampling 30j, milestoning [31], and Markov state 



models 



32l | have been developed in the last decade. 



In this paper, we pursuit another approach by discretizing the OM action 

N-l 



S 



AU 



C 



where n = 1/Ai, X\ = x-^xn 



i=l 

Xfi n and 



+ Yl o( X *+l -Xi) 2 + W(Xi 



(3) 



W(x) = i- 
Ik 



( dx J 



2Dd 2 U 



(4) 



The statistical weight is P(xi,X2, ■ ■ ■ ,xn) oc e~ s ^ 2D , and the path search problem is there- 
fore mapped onto the statistical mechanics of a polymer under the effective potential function 
W(x) 28 1 . (This isomorphology is similar to that between a quantum particle and a ring 



34fl but has not been 



polymer 33].) This strategy was proposed in the original TPS paper 
completely worked out. The previous studies using this type of action-based method have 
mostly focused on a single most probable path during conformational changes of biomolecules 



35H38] . Here, we try to complete the previous studies and provide a path ensemble at finite 



temperatures using the action-based method, which should be relevant for slow and diffusive 
biomolecular reactions. 

To generate a path ensemble, we need to assume an initial path, from which the path 
sampling simulation starts. If there is a barrier in path space, however, the generated 
path cannot move to the most probable path because "path trapping" occurs in the basin 
around the initial path (as in the case of "configuration trapping" for the folding problem 
of proteins). Even worse, for complex systems with rugged energy surfaces, the path to be 
calculated may not be unique, and multiple paths can coexist. This is the situation we want 
to address in this paper. In the previous study using the nudged-elastic band, simulated 



annealing was employed to avoid this "path trapping" problem [18] . This is a nice strategy 
to sample multiple minimum energy paths, however, a path ensemble at a finite temperature 
cannot be obtained. To generate a canonical ensemble for a path, we propose to introduce 



39 



4Q|, 



one of the generalized ensemble methods, the replica exchange method (REM) 
in path space. The combination of REM with transition path sampling has been already 
proposed 4lM43] for rare and fast chemical reactions, however, our aim is to sample rare, 
slow, and diffusive processes, and the application of the OM action formalism should be 
more appropriate. 



This paper is organized as follows. In Sec. 011 we derive an overdamped Langevin equation 
for a path such that the path ensemble is generated according to the weight determined by 
the OM action. In parallel to the technique employed in the path integral simulations of 



quantum systems, we use Fourier components of a path as "dynamic" variables |44j. We 
then combine this method with REM to improve the sampling efficiency in path space. 
In Sec. IHIt we test the numerical performance of our new method. To this end, a two- 



dimensional model potential due to Bolhuis 



43], involving bifurcated reaction pathways, is 



numerically examined in detail. The results are finally compared with those of TPS and the 
equilibrium theory. In Sec. IIV} we summarize the paper and discuss the connection between 
our method and other related ones and future development for biomolecules. 



II. THEORY 



A. Onsager-Machlup action for multi-dimensional systems 



To describe slow and diffusive processes, we start from an overdamped Langevin equation 
for an M-dimensional system using mass-weighted Cartesian coordinates (xi, ■ ■ ■ ,xm)'- 

1 dU 



X r 



+ y/2D aVa (t) 



(5) 



J a dx a 

where 7 a is an intrinsic friction coefficient, i] a (t) is a Gaussian-white noise satisfying 
(r]a(t)ria'(t')) = 8 aa i5(t - t'), and 

A, = ^ (6) 

is imposed by the fluctuation-dissipation theorem. The corresponding Fokker-Planck equa- 
tion [45] is 

dP({x a >},t) d (\ dU „ /f .A d 2 n/f . , 



Of 



dx a \7q dx c 



The path-integral representation of the propagator (Green's function) is written as 29l. |38| 

P(x fin |x ini ;t)oc / ,m D I ( s)e -«[^/^ ) (8) 

J x(0)=Xi n [ 

where H is the OM action defined by the multidimensional extension of Eq. ([2]), multiplied 
by 7 a /2 and summed over a. Its discretized form is written as 



^ AC/ K 

i=i 



M 



a=l 



^ - X a)i ) 2 + V e g({x ayi }) 



(9) 



where AU = U{x^ D ) — U(x- U) i) and the effective potential is defined as 



V eS ({x a }) 



with 



M 

E 

a=l L 



u 2 

2 x o 



k B T, 



2 Oi&OL 



To 



2At 



(10) 



(11) 



being an effective frequency determined by the friction and the time interval for the path 
discretization At. (Hereafter U X ,U XX represent the first and second derivatives of U with 
respect to x.) Note that the path ensemble generated by the above OM "Hamiltonian" is 
isomorphic to the canonical ensemble of an M x (N — 2) dimensional (polymer) system. 



B. Path sampling in Fourier space 

For simplicity, we first consider a one-dimensional system. According to Cho, Doll, and 
Freeman 44j, the Fourier transform of a path connecting x\ and x^ is defined as 

V k=2 ^ ' k=2 

where qk is the Fourier component of the path x i} and is £1 reference path, which may be 
an initially guessed path. To generate a path ensemble with the weight determined by the 
OM Hamiltonian T-t, Eq. ()9]), the following overdamped Langevin dynamics can be employed 

qk = -L k ^(/3 T H) + ^/2T kVk (t), (13) 
OQk 

with 

(Vk(tMt')) = 8 kl S(t-t'), (14) 

where L k is an Onsager coefficient representing friction, and 1/HbPt is the temperature 
of a thermostat, not necessarily equal to the system temperature T. Using the corre- 
sponding Fokker-Planck equation, the stationary point is shown to be the equilibrium state 
P eq ({q k }) oc exp(— f3 T l-L) in (discretized) path space. Note that this Langevin dynamics is 
different from the original Langevin dynamics, Eq. (jSJ), which was motivated from a physical 
consideration. 
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Substituting the particular form of the OM Hamiltonian, Eq. fl9]), into the overdamped 
Langevin equation, we have 

N-l 

4k = 



i=2 l 
— = -W + Xi-! -2Xi) + 



dxi 



dV eS (x) 



l -uu - k -^u 



cte .U-^ " 4a; 2 ^' 

where a; is the one-dimensional analog of Eq. (II ip . Finally we obtain 



i=2 



where 



N-l 



G k = W 2 /3 T L fc ^ Mjfc (x l ( ° ) 1 + xfJ 1 -2xf ) ), 

n(k - iy 



i=2 



2ui firLk 1 — cos 



N-l 



By taking 



1 



1 — cos 



n(k - 1) 



(15) 

(16) 
(17) 

(18) 

(19) 
(20) 

(21) 



\At/3 T uj 2 \~ ~~~ N-l 
the rate becomes constant: = To = 2/(XAt) where A is a dimensionless parameter to 

adjust a timescale for relaxation. This particular form of the Onsager coefficient was chosen 

because the time step to solve the above series (k — 2, 3, ■ • • , N — 1) of Langevin equations 

may be determined by a single timescale 1/T . 

The Fourier-path Langevin dynamics for the M-dimensional system becomes 



N-l 



i=2 



dV eS ({x a }) 



dx r 



+ y^T k r] a:k (t), (22) 



where 



dV cS ({x a }) 



M 

£ 

)3=1 



knT^ 



U Tj —U 



H 1 



(23) 



and L k is chosen as in Eq. (|2~T|) with w = w Q . Practically, the hessian and its derivative 
above are the bottleneck of this computation, hampering its application to large molecular 
systems. Recently, however, the symmetric OM action, which only uses the force to calculate 
the OM action, has been devised by Miller and Predescu [46]. Such a "low cost" action can 
be utilized for the future application of this method to real molecular systems. 



C. Replica exchange in path space 



In the previous subsection, we have obtained a primitive way to generate a path ensemble 
by solving the overdamped Langevin equation for Fourier components of a path. However, 
path sampling efficiency becomes an issue in biomolecular simulations due to the extremely 
rugged-energy surface. In such a case, the sampled path often gets trapped around the 
initially assumed path, leading to unphysical results and incorrect predictions on the reaction 
mechanism. This situation is in parallel with the configuration sampling problem of complex 
molecules such as peptides or proteins in a rugged-energy landscape. 

Genera 



change 
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ized ensemble methods such as multi-canonical sampling [47], |48j] or replica ex- 
] (parallel tempering) are well-known to solve this problem. We here combine 
one of these methods, the replica exchange method (REM), with the OM method because 
of its simplicity for implementation. It is easy to derive the following acceptance probability 



for exchange of two replicas 

P acc = min(l,exp{A}), (24) 
A = (0-p')(H-H'), (25) 

where (/3,1-L) and (^',H') are the inverse temperature and the OM Hamiltonian for each 
replica. The difference with the conventional REM 0, is that each replica is associated 
with a path (see Fig. [I]), and that Eq. (Q is employed for the effective Hamiltonian H. As 
usually done, we prepare several replicas for a path and, during the Langevin dynamics in 
our case, exchange two "neighboring" paths at the same time with the above probability 



40j. If the number of replicas and the exchange rate between replicas are both sufficiently 
large to sample the whole "path" space, an appropriate path ensemble with temperature T 
should be obtained by collecting path data indexed by the same T. 
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FIG. 1: Schematic drawing of the proposed replica-exchange method for path sampling using the 
Onsager-Machlup action. Each column represents the OM dynamics of a single path. Different 
path replicas with different inverse temperatures /3j (i = 1,2, •• • ,n) are exchanged according to 
the Metropolis criterion, Eqs. (|24p and (|25p . 

III. RESULTS 

A. Model and methods 




We test our method by using Bolhuis' two-dimensional potential [431] 



U(x,y) = -3e-°- 25 ^ 2 -y 2 -3e-°- 25 ^ 2 -y 2 

+ JL( .0625a; 4 + y 4 ) + SeT " 0081 * 4 - 4 *' 2 
1800 V y 1 

+2e -1.5( X -b)^(y-lf + 2ae _i.5(* + 6) 3 -(y+l) a _ (26 ) 

This potential has two minima at (±4.3, 0.0) which are interconnected by two reaction 
pathways via two saddle points (0.0, ±2.3). These pathways are separated by a large energy 
barrier centered at (0.0,0.0). The parameters are taken as a = 1 and 6 = 0, resulting in 
a symmetric shape of the potential (see Fig. [5]). This potential was chosen because it is 
a minimal model with multiple pathways, representing the local character of the rugged 
energy landscape of biomolecules. Our aim is to sample an ensemble of paths connecting 



10 



6-4 



x 



FIG. 2: Left: Bolhuis' potential defined by Eq. (|26[) with a = 1 and 6 = (symmetric case). 
There are two stable regions (basins) around (x, y) = (±4.3, 0) and a large barrier separating two 
pathways connecting the two basins. Right: Contour plot of the Bolhuis potential. The red and 
green vertical lines represent the exit positions used later in the analysis. The blue curve connecting 
two minima represents a minimum energy path and the intersections with the red and green lines 
are denoted as triangles. 

two minima. 

We need to determine the parameters for the overdamped Langevin equation, Eq. (JSJ). 
We first take 7 Q = 1, i.e., the characteristic timescale of this Langevin dynamics is unity. 
From the Einstein relation, Eq. ([6]), the diffusion constant is equal to the thermal energy, 
ksT. Here we set ksT = 1.25, since this is an interesting case where barrier crossing can 
easily occur in the x direction, but not in the y direction. 

Next we describe how to discretize a path in the Bolhuis potential. The discretization of 
a path can be characterized by the time interval At and the number of the discretized path 
(the number of "beads") N, i.e., the total time of the path is t to t = NAt. To estimate the 
reasonable value for t tot , we consider free diffusion between the two minima with length L, 
where the diffusion timescale is t^is = L 2 /2D. In our case, L ~ 2 x 4.3 = 8.6 and D = 1.25, 
and thus t^m — 30. However, if the potential surface is concave in the transition state region, 
the motion can be quicker than free diffusion, and the actual time for barrier crossing can 
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be shorter than this estimate. In this paper, we chose t to t = 12. The next question is how 
to choose At or iV while NAt = 12. After several trial and error, we found that At = 0.05 
and N = 240 are appropriate for this Bolhuis potential. 

We further need to determine the time step for the Fourier-path Langevin dynamics, 



Eq. We used the simplest Euler scheme to solve this dynamics [49|, and found that 

A = 10 5 and AtF = 0.5 are reasonable. Note that the Langevin dynamics for a path is carried 
out just for sampling and the timescale determined by these parameters have no physical 
meaning. For comparison, we also numerically solve the original overdamped Langevin 
equation, Eq. ([S]), with the same Euler scheme. We have confirmed that Ati = 0.003 is 
sufficient to recover the potential function through the formula U(x, y) = —k B T log P(x, y) + 
C, where P(x,y) is the distribution function calculated by the Langevin trajectory and C 
is a constant. 

Finally we describe how to choose the parameters for REM. The most important prerequi- 
site for REM is that the distribution of H should have sufficient overlap between neighboring 



replicas [39 



40| . which is written as 

\{n)i+i - (U)i\ < (AH) i+ x + (AH)i. (27) 

A rough estimate may be obtained from an independent particle approximation and equipar- 
tition principle. For the spring contribution in Eq. (jHJ), the number of particles is approxi- 
mately N x M, and the thermal energy k B T/2 is given to each degree of freedom, the partial 
average value of the Hamiltonian % is therefore N Mk B T '/2. Actually there is a contribution 
from the effective potential term V e g, which might be also approximated as an A^-tiple con- 
nected harmonic spring, then the average of the total Hamiltonian may be (H) = NMk B T. 
Assuming that they are all uncorrelated Gaussian random numbers, the fluctuations of the 
Hamiltonian becomes (AH) = y/2N MksT '. Defining / = Tj + i/Tj > 1, the criterion for the 
neighboring distributions to overlap is written as 

NM(f - 1) < V2NM(f + 1), (28) 

and we find 

1< / < 1 + ; 2 . (29) 

yfNM/2-l 

In our case of the Bolhuis potential, NM ~ 480, the right-hand side thus becomes 1.14. In 
this paper, we prepare eight replicas by taking / = 1.14, hence the corresponding temper- 
atures are k B T = 3.60, 3.10, 2.66, 2.29, 1.97, 1.69, 1.46, 1.25. We have checked that the eight 
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FIG. 3: Path ensemble generated by the Fourier-path dynamics without REM (left) and with REM 
(right). The diffusion constant or temperature is D = ksT = 1.25. 200 paths are superimposed 
for each panel with the interval of 30 x Atp- The initial guess is also drawn as a green line. 

replicas are mixed up well and the equilibrium path distributions are sampled appropriately 
at the respective temperatures. 



B. Numerical results 



In Fig. |3l we show the path ensembles generated by the Fourier-path dynamics with 
(right) or without (left) the use of REM. An initially guessed path (green line in Fig. [3]) 
was shifted from a linearly interpolated path to avoid the energy being too high. Starting 
from the initial path, the generated path gradually spreads out in path space. However, in 
this "timescale" (200 x 30 x Atp), the path sampling without REM is not sufficient because 
the distribution turns out not to be symmetric (see the left of Fig. [3]). This indicates that 
the simple Fourier-path dynamics does not converge fast enough when there are multiple 
pathways: this is the situation of "path trapping" that we may encounter in path sampling 
simulations of biomolecules. For such a case, REM should be employed to improve the 
sampling in path space. As shown in the right of Fig. [3j the Fourier-path dynamics with 
REM can explore larger path space, resulting in a more symmetrical distribution as expected. 

From the path ensemble thus generated, we can obtain both nonequilibrium (such as 
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kinetic rates) and equilibrium properties (such as free energy profiles). Since we want to 
examine the nonequilibrium properties of the path ensemble, we employ an exit distribution 



50( as a measure of such nonequilibrium properties. Note that the exit distribution is 
conceptually different from a free energy surface which is defined for equilibrium systems. 
For our model system, the exit distribution is defined as a histogram of accumulated y when 
a trajectory first hits a cross section x = x ex it- In Fig. HJ we show the numerical result of the 
exit distributions for two exit positions: (a) one is on the transition state (TS): a; ex it = 0, 
corresponding to the red line in Fig. El and (b) the other is apart from the TS: x exit = —1.33, 
corresponding to the green line in Fig. EJ Because it is located on TS, case (a) is expected to 
show stronger nonequilibrium behavior than case (b). To check the correctness of our result, 
we also show the exit distributions calculated by transition path sampling (TPS) using the 
algorithm due to Crooks and Chandler 50] . For the TPS calculations, the total time was 
t to t — 12 as we used in the Fourier-path dynamics calculation, and a small cutoff distance 
(0.05) was introduced to define the reactant and product states around two minima. We see 
that two exit distributions calculated by the two methods (at the different exit positions) 
agree well within the statistical error. This result indicates that our method using the OM 
action with REM can be alternative to TPS. From computational points of view, for this type 
of "small" system, TPS works best because it can easily sample the whole path space using 



the Croo 



ts-Chandler algorithm. The algorithm due to Eastman, Gronbech- Jensen, and 



Doniach 37J is almost comparable to ours because both methods use the same information 
of the system such as the derivative of a hessian matrix. Note, however, that our intention 
is to present an efficient path sampling algorithm for large systems. This approach of path 
sampling is considered to be more feasible for slow processes in large systems. 

Next, the exit distributions are compared to the equilibrium probability function (PDF) 
at the exit points, which is calculated as 

P eq {y) oc exp{-(3U(x exit ,y)} (30) 

after normalization along the y direction. It is interesting to see that although the PDF is 
conceptually different from the exit distribution, they look similar to each other at x cx it = 
— 1.33 (case (b), see the right of Fig. HJ). On the other hand, the exit distribution is quite 
different from the PDF at x ex it = 0.0 (case (a), see the left of Fig. HJ), which should be 
attributed to nonequilibrium effects. We have also computed the exit distributions using 
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FIG. 4: Exit distributions on the transition state (x ex n = 0, left, corresponding to the red line on 
the right panel of Fig. [2|) and apart from the TS (x ex it = —1-33, right, corresponding to the green 
line on the right panel of Fig. [2]) calculated by the OM method with REM (red) and TPS (green) 
with D = ksT = 1.25. The equilibrium probability distribution function (PDF) is also shown 
(blue). 

direct Brownian dynamics starting from one minimum (data not shown). Interestingly, the 
result also gives a similar exit distribution for case (a) but not for case (b). This indicates 
that the exit distribution does depend on the final state as well as the initial state, and for 
case (b) such an effect can be significant. As such, our method or TPS should be carefully 
compared with direct Brownian dynamics. 

The calculations done have been so far successful, and we now discuss practical aspects 
of the path sampling simulations. In our approach, one of the important parameters that 
affect the computational effort is the length for path discretization At. It is thus useful to 
mention how and why the calculation fails when we take larger At. In the left of Fig. [51 we 
show the result of the exit distribution with At = 0.1 and N = 120 for ) while fixing; 

the total time t tot = 12 (the original setting was At = 0.05 and N = 240). Compared with 
the left of Fig. HJ one can see that spurious peaks appear around x = ±1.0. 

To understand these peak structures, we carefully examine the effective potential, 
Eq. ( TTU1) . This potential can be decomposed into two parts: one is the zero-temperature 
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FIG. 5: Left: Exit distributions calculated by our method and TPS at the exit position x ex it = 
with ksT = 1.25. For the OM calculation, At = 0.1 and N = 120 were chosen for the path 
disretization. The probability distribution function at the exit position is also shown as a reference. 
Right: Cross sections of the effective potential (red) and its zero-temperature component (green) 
at the exit position x GX ;t = with fcgT = 1.25. The cross section of the potential energy function 
is also shown in blue. 

component v}g({x a }) and the other the finite-temperature component v}^ ({x a }), that is, 

V eS (M) = V<£\{x a }) + Vg{{x a }), (31) 

M 

V£\{x a }), = £tt-^L (32) 

a=l a 
M 

vS\{x a }) = -k B Tj2^U XaXa . (33) 

a=l a 

It is easily seen that the finite-temperature component becomes zero when T = 0. More 
importantly, the finite-temperature component prefers a path with positive and large cur- 
vatures. This is because an ensemble of trajectories becomes stable in the region where a 
curvature is large. (If the curvature is negative, the ensemble is unstable and spreading.) 

In the right of Fig. [5j we show the cross sections at the exit points for and V e g. From 
this figure, we can clearly see that the zero-temperature component does not have a sufficient 
ability to drag the path around y ~ ±1.0. We thus conclude that the peak structures in 
the exit distribution stem from the finite-temperature component of the effective potential. 
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Hence we need to be careful when we deal with a path ensemble with large At because some 
spurious features such as the peak structures might appear. 



IV. CONCLUDING REMARKS 

As a novel approach to sample diffusive paths based on the Onsager-Machlup action, 
we have proposed Fourier-path Langevin dynamics. To achieve powerful sampling in path 
space and to avoid the problem of "path trapping" around an initially guessed path, we 
have also suggested to combine this scheme with a powerful sampling technique, the replica 
exchange method. Using the two dimensional model potential due to Bolhuis, the validity of 
our method has been confirmed by the numerical comparison to the conventional transition 
path sampling. We have also identified an interesting nonequilibrium path ensemble near 
the transition state of the model system, which is different from an equilibrium distribution. 

We are now in a position to discuss the relation of our method with other path search 
or path sampling methods. First, our main concern is the effects of temperature on a p ath, 



so the path search methods at a finite temperature such as MaxFlux methods 5l|, |52| . 



;he temperature-dependent reaction coordinate 



54 



531 ]. or finite-temperature string methods 



55] have strong similarity with our method. The difference is that a path generated 
by our method still holds nonequilibrium properties of the path and there is no assumption 
such as the existence of local equilibrium. It is interesting and important to investigate how 
the path ensembles calculated by different methods are actually different. 

In this paper, for simplicity, we assumed overdamped Langevin dynamics, but this restric- 
tion can be easily relaxed. We can employ a modified action derived from the underdamped 
Langevin dynamics instead of the OM action [56] . This strategy was successfully used for 
dynamic reweighting of a trajectory for a model system [53]. When we apply this method 
to real molecular systems, we need to judge which Langevin dynamics (overdamped or un- 
derdamped) is more appropriate. 

Of course, Eq. (fl3|) is not the only dynamics that can generate a canonical ensemble of a 
oath, e.g., one may use NVT molecular dynamics with Nose-Hoover thermostat techniques 
|59j]. In principle, a much simpler algorithm can be used such as a Monte Carlo (MC) 



move in path space. That is, we devise a certain MC move x — > x', and accept or reject the 



move using the Metropolis criterion 42|. However, we need a clever move when we apply 
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the MC algorithm to biomolecular systems, otherwise the move is rarely accepted. The 
reason for this deficiency is the same as why a simple MC move does not work for sampling 
in biomolecular configuration space, i.e., a MC move without consideration of a molecular 



configuration causes a high energy state which will be rejected [60j. If we can devise a clever 
move in path space, it would be a powerful alternative to sample path space. Combining 
the OM method with the other generalized ensemble methods such as the multi-canonical 



method 



47 



48j or the Tsallis ensemble method 6l|, |62| is also possible and promising. 
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